<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Minimize stopband ripple of a linear phase lowpass FIR filter</title>
<link rel="canonical" href="http://cvxr.com/cvx/examples/filter_design/html/fir_lin_phase_lowpass_min_ripple.html">
<link rel="stylesheet" href="../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Minimize stopband ripple of a linear phase lowpass FIR filter</h1>
Jump to:&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#source">Source code</a>&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#output">Text output</a>
&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#plots">Plots</a>
&nbsp;&nbsp;&nbsp;&nbsp;<a href="../../index.html">Library index</a>
</div>
<div id="content">
<a id="source"></a>
<pre class="codeinput">
<span class="comment">% "Filter design" lecture notes (EE364) by S. Boyd</span>
<span class="comment">% (figures are generated)</span>
<span class="comment">%</span>
<span class="comment">% Designs a linear phase FIR lowpass filter such that it:</span>
<span class="comment">% - minimizes the maximum passband ripple</span>
<span class="comment">% - has a constraint on the maximum stopband attenuation</span>
<span class="comment">%</span>
<span class="comment">% This is a convex problem.</span>
<span class="comment">%</span>
<span class="comment">%   minimize   delta</span>
<span class="comment">%       s.t.   1/delta &lt;= H(w) &lt;= delta     for w in the passband</span>
<span class="comment">%              |H(w)| &lt;= atten_level        for w in the stopband</span>
<span class="comment">%</span>
<span class="comment">% where H is the frequency response function and variables are</span>
<span class="comment">% delta and h (the filter impulse response).</span>
<span class="comment">%</span>
<span class="comment">% Written for CVX by Almir Mutapcic 02/02/06</span>

<span class="comment">%********************************************************************</span>
<span class="comment">% user's filter specifications</span>
<span class="comment">%********************************************************************</span>
<span class="comment">% filter order is 2n+1 (symmetric around the half-point)</span>
n = 10;

wpass = 0.12*pi;        <span class="comment">% passband cutoff freq (in radians)</span>
wstop = 0.24*pi;        <span class="comment">% stopband start freq (in radians)</span>
atten_level = -30;      <span class="comment">% stopband attenuation level in dB</span>

<span class="comment">%********************************************************************</span>
<span class="comment">% create optimization parameters</span>
<span class="comment">%********************************************************************</span>
N = 30*n+1;                            <span class="comment">% freq samples (rule-of-thumb)</span>
w = linspace(0,pi,N);
A = [ones(N,1) 2*cos(kron(w',[1:n]))]; <span class="comment">% matrix of cosines</span>

<span class="comment">% passband 0 &lt;= w &lt;= w_pass</span>
ind = find((0 &lt;= w) &amp; (w &lt;= wpass));   <span class="comment">% passband</span>
Ap  = A(ind,:);

<span class="comment">% transition band is not constrained (w_pass &lt;= w &lt;= w_stop)</span>

<span class="comment">% stopband (w_stop &lt;= w)</span>
ind = find((wstop &lt;= w) &amp; (w &lt;= pi));  <span class="comment">% stopband</span>
Us  = 10^(atten_level/20)*ones(length(ind),1);
As  = A(ind,:);

<span class="comment">%********************************************************************</span>
<span class="comment">% optimization</span>
<span class="comment">%********************************************************************</span>
<span class="comment">% formulate and solve the linear-phase lowpass filter design</span>
cvx_begin
  variable <span class="string">delta</span>
  variable <span class="string">h(n+1,1)</span>;

  minimize( delta )
  subject <span class="string">to</span>
    <span class="comment">% passband bounds</span>
    Ap*h &lt;= delta;
    inv_pos(Ap*h) &lt;= delta;

    <span class="comment">% stopband bounds</span>
    abs( As*h ) &lt;= Us;
cvx_end

<span class="comment">% check if problem was successfully solved</span>
disp([<span class="string">'Problem is '</span> cvx_status])
<span class="keyword">if</span> ~strfind(cvx_status,<span class="string">'Solved'</span>)
  <span class="keyword">return</span>
<span class="keyword">else</span>
  <span class="comment">% construct the full impulse response</span>
  h = [flipud(h(2:end)); h];
  fprintf(1,<span class="string">'The optimal minimum passband ripple is %4.3f dB.\n\n'</span>,<span class="keyword">...</span>
            20*log10(delta));
<span class="keyword">end</span>

<span class="comment">%********************************************************************</span>
<span class="comment">% plots</span>
<span class="comment">%********************************************************************</span>
figure(1)
<span class="comment">% FIR impulse response</span>
plot([0:2*n],h',<span class="string">'o'</span>,[0:2*n],h',<span class="string">'b:'</span>)
xlabel(<span class="string">'t'</span>), ylabel(<span class="string">'h(t)'</span>)

figure(2)
<span class="comment">% frequency response</span>
H = exp(-j*kron(w',[0:2*n]))*h;
<span class="comment">% magnitude</span>
subplot(2,1,1)
plot(w,20*log10(abs(H)),[wstop pi],[atten_level atten_level],<span class="string">'r--'</span>);
axis([0,pi,-40,10])
xlabel(<span class="string">'w'</span>), ylabel(<span class="string">'mag H(w) in dB'</span>)
<span class="comment">% phase</span>
subplot(2,1,2)
plot(w,angle(H))
axis([0,pi,-pi,pi])
xlabel(<span class="string">'w'</span>), ylabel(<span class="string">'phase H(w)'</span>)
</pre>
<a id="output"></a>
<pre class="codeoutput">
 
Calling sedumi: 884 variables, 606 equality constraints
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
Split 12 free variables
eqs m = 606, order n = 860, dim = 934, blocks = 38
nnz(A) = 1330 + 6814, nnz(ADA) = 1286, nnz(L) = 946
Handling 24 + 0 dense columns.
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            4.89E-02 0.000
  1 :   6.36E+00 9.87E-03 0.000 0.2019 0.9000 0.9000  -0.43  1  1  2.9E+00
  2 :   1.94E+00 5.09E-03 0.000 0.5157 0.9000 0.9000   4.74  1  1  4.2E-01
  3 :   1.15E+00 3.34E-03 0.000 0.6558 0.9000 0.9000   9.02  1  1  6.5E-02
  4 :   1.08E+00 1.96E-03 0.000 0.5865 0.9000 0.9000   2.40  1  1  2.8E-02
  5 :   1.06E+00 1.12E-03 0.000 0.5714 0.9000 0.9000   1.65  1  1  1.4E-02
  6 :   1.05E+00 4.73E-04 0.000 0.4224 0.9000 0.9000   1.41  1  1  5.2E-03
  7 :   1.05E+00 1.18E-04 0.000 0.2501 0.9000 0.9000   1.19  1  1  1.2E-03
  8 :   1.05E+00 4.38E-05 0.000 0.3703 0.9000 0.9000   1.05  1  1  4.5E-04
  9 :   1.05E+00 7.92E-06 0.000 0.1810 0.9000 0.0000   1.01  1  1  2.6E-04
 10 :   1.05E+00 1.80E-06 0.000 0.2274 0.9166 0.9000   1.01  1  1  7.3E-05
 11 :   1.05E+00 5.21E-07 0.000 0.2895 0.9000 0.8557   1.00  1  1  2.0E-05
 12 :   1.05E+00 1.98E-07 0.000 0.3797 0.9035 0.9000   1.00  2  2  7.9E-06
 13 :   1.05E+00 6.25E-08 0.000 0.3158 0.9000 0.8355   1.00  2  3  2.4E-06
 14 :   1.05E+00 1.63E-08 0.000 0.2609 0.9000 0.9000   1.00  4  4  6.4E-07
 15 :   1.05E+00 3.20E-09 0.000 0.1963 0.9000 0.9000   1.00  5  5  1.2E-07
 16 :   1.05E+00 2.20E-10 0.000 0.0687 0.9900 0.9900   1.00  7  7  8.6E-09

iter seconds digits       c*x               b*y
 16      0.3   Inf  1.0515780150e+00  1.0515780156e+00
|Ax-b| =   4.8e-08, [Ay-c]_+ =   3.5E-09, |x|=  1.2e+01, |y|=  1.2e+00

Detailed timing (sec)
   Pre          IPM          Post
1.000E-02    2.700E-01    0.000E+00    
Max-norms: ||b||=1, ||c|| = 1,
Cholesky |add|=6, |skip| = 2, ||L.L|| = 1.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +1.05158
Problem is Solved
The optimal minimum passband ripple is 0.437 dB.

</pre>
<a id="plots"></a>
<div id="plotoutput">
<img src="fir_lin_phase_lowpass_min_ripple__01.png" alt=""> <img src="fir_lin_phase_lowpass_min_ripple__02.png" alt=""> 
</div>
</div>
</body>
</html>